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In this work the free energy of solid phases is computed for the Lennard- Jones potential and for a 
model of NaCl. The free energy is evaluated through the Einstein crystal/molecule methodologies 
using the Molecular Dynamics programs: GROMACS and LAMMPS. The obtained results are 
compared with the results obtained from Monte Carlo. Good agreement between the different 
programs and methodologies was found. The procedure to perform the free energy calculations for 
the solid phase in the Molecular Dynamic programs is described. Since these programs allow to 
study any continuous intermolecular potential (when given in a tabulated form) this work shows 
that for isotropic potentials (describing for instance atomic solids or colloidal particles) free energy 
calculations can be performed on a routinely basis using GROMACS and/or LAMMPS. 



In 1984 Frenkel and Ladd proposed the "Einstein crys- 
tal method" , a novel scheme to compute the free energy 
of a solid flj: the method is based on a thermodynamic 
integration of the Helmholtz free energy in the canonical 
ensemble along a reversible path between the system of 
interest and an ideal Einstein crystal with the same struc- 
ture as the real solid, whose Helmholtz free energy can be 
analytically computed. In the ideal Einstein crystal, par- 
ticles are attached to their lattice positions via harmonic 
springs (of constant As). More recently, some of us have 
proposed the "Einstein molecule method" 0, Q , a vari- 
ant of the Einstein crystal to compute the free energy 
of molecular solids. The main difference between both 
methods is the choice of the reference system. In the 
Einstein crystal the reference system is an ideal Einstein 
crystal with the constraint of the center-of-mass of the 
system kept fixed (to avoid a quasi-divergence of the in- 
tegral of the free energy change from the reference system 
to the real solid). In the Einstein molecule the reference 
system is an ideal Einstein crystal with the constraint of 
the position of one particle kept fixed. 

The Helmholtz free energy A so i (T, V) computed with 
the Einstein crystal/molecule calculations can be writ- 
ten as A sol (T, V) — A (T, V) + AA± (T, V) + AA 2 (T, V) , 
where Ao is the free energy of the reference system (in- 
cluding corrections for the fixed point), whose analytical 
expression is slightly different in the Einstein crystal and 
Einstein molecule, AAi is the free energy difference be- 
tween the ideal Einstein crystal and the Einstein crystal 
in which particles interact through the Hamiltonian of 
the real solid ("interacting" Einstein crystal) and AA2 
is the free energy difference between the interacting Ein- 
stein crystal and the real solid [4j. The expression for 
both AAi and AA2 is the same in both methods although 
its value may be different since the fixed point (center- 
of-mass or a reference particle) is different (a detailed 
description of these terms is provided as Supplementary 
Material[|). In any case since the free energy of a solid 
is uniquely defined, its value should not depend on the 
methodology used to compute it 0, ||J . 

Besides the Einstein crystal methods other methods 



have been developed in the last decade to estimate the 
free energy of solids (or to determine the fluid-solid equi- 
libria), such as the Phase Switch method @ and an inte- 
gration path to calculate the Gibbs free energy differ- 
ence between any arbitrary solid and liquid proposed 
by Grochola 0]. In spite of this progress the number 
of groups performing free energy calculations for solids 
is still small Solid free energy calculations can 

be easily adapted into Monte Carlo scheme, but it re- 
quires one to write a la carte MC codes. Therefore it 
seems of interest to consider if these calculations can 
be implemented in widely used open source Molecular 
Dynamics simulation programs such as GROMACS [ll[ 
or LAMMPS [12j. As far as we are aware, free en- 
ergy calculations of solids are not yet explicitly imple- 
mented in such packages and only recently free energy 
calculations for alloys of copper and zirconium [l3j us- 
ing LAMMPS through the Einstein crystal method have 
been reported. Since GROMACS and/or LAMMPS in- 
corporate the possibility of using harmonic potentials re- 
straining the position of atoms (or even freezing the po- 
sition of atoms) and incorporate thermostats that can 
correctly treat harmonic oscillators (as for instance the 
new thermostat of Bussi et al.[14|), it seems possible to 
perform Einstein crystal/molecule calculations for atomic 
solids using these programs. In this note we shall use 
GROMACS and LAMMPS to compute the free energy 
of solids of particles interacting via isotropic potentials 
using either the Einstein crystal or the Einstein molecule 
method. In particular we will compute the free energy of 
Lennard- Jones (LJ) and of Sodium Chloride modeled via 
the Joung-Cheatham(JC)-NaCl potential [l5[ (using the 
SPC/E set) which describes ion- ion interactions through 
a LJ potential plus a Coulomb term. To confirm the re- 
sults obtained with GROMACS (G) and LAMMPS (L), 
calculations will be also performed using a Monte Carlo 
code (MC). 

We have carried out free energy calculations for a 
256 LJ-Argon spherically truncated and shifted (STS) 
at r c = 2.7a and for a 1372 LJ spherically truncated 
(ST) at r c = 5a. Results of the free energy calculations 
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TABLE I. Free-energies of solids as obtained from Ein- 
stein Crystal (EC) and Einstein Molecule (EM) method- 
ologies using MC or MD (G;L). N is the number of parti- 
cles, r c is the cut-off distance of the LJ contribution. Re- 
sults for the LJ/STS and LJ/ST systems were obtained at 
T* = k B T/e = 2 and p* = pa s = 1.28 for the fee structure. 
Results for the NaCl solid were obtained for the JC-NaCl 
model at T = 298 K and a volume of V = 24.13 nm 3 . Free 
energies are in NIcbT units. 



are presented in Table Q] For the LJ/STS system A , A 
Ai, and A A 2 obtained from MC and from MD (GRO- 
MACS/LAMMPS) are in very good agreement, with a 
free energy difference lower than 0.03 NksT (typical un- 
certainties in calculations of free energy of solids being 
of 0.05 NksT). The choice of a LJ spherically truncated 
and shifted (STS) avoids the subtle issues arising when 
comparing results obtained by MC and MD for a spher- 
ically truncated (ST) potential!!, [H, [T^. In the LJ/ST 
truncated at r c — 5<r, the free energy results obtained 
with MC and MD agree quite well with each other and 
with previous calculations for the same system size and 
thermodynamic state @, [H, [lj| . It is clear that the mag- 
nitude of problems arising with the discontinuity of the 
potential at the cut-off are also quite small when the cut- 



off distance is large enough (« 5a). 

Notice that although for a certain system size and 
Hamiltonian the computed free energy of a solid is 
unique, in general the free energy of solids changes with 
the system size. As a suggestion, in order to have a rea- 
sonable estimate of the free energy of a solid, we recom- 
mend to use a relatively large system size (with more than 
1000 molecules) and to average properties over around 
10 4 — 10 5 independent configurations for evaluating AAi 
and runs of about 10 ns for AA2. Moreover, the use of 
a large cut-off (above 4.5tr) is recommended. In all LJ 
systems, both Einstein molecule and Einstein crystal give 
the same value of the free energy A so i , although the single 
contributions are slightly different in both methodologies 
(as they should be). In the Supplementary Material, we 
provide further details on the implementation |5j, needed 
for GROMACS and LAMMPS to easily compute A , A 
Ai, and A A 2 . 

To show that the methodology also works for more 
complex systems where particles interact via an isotropic 
potential, we calculate the free energy of NaCl using the 
Joung- Cheatham potential [Hj]. Results are presented in 
Table H We simulate a 1000 ions NaCl solid (at 298 K 
and the equilibrium density of the model at 1 bar). The 
free energies computed from MC and GROMACS agree 
reasonably well (the difference been of 0.03 NknT) and 
are in good agreement with previous calculations [20] . 

To conclude, we have demonstrated that it is possible 
to compute the free energy of atomic solids using GRO- 
MACS and LAMMPS for systems interacting with spher- 
ical potentials such as Lennard- Jones, Yukawa, Morse 
and any continuous potential since both packages incor- 
porate the possibility of using a tabulated numerical in- 
termolecular potentials. We do hope that this work stim- 
ulate more groups to perform free energy calculations for 
solids. For future work it would be useful to analyze if 
for molecular fluids the Einstein crystal/molecule calcu- 
lations could also be performed with GROMACS and/or 
LAMPPS. 
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